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ABSTRACT 

In a semi-numerical model of reionization, the evolution of ionization fraction is simulated 
approximately by the ionizing photon to baryon ratio criterion. In this paper we incorporate a 
semi-analytical model of galaxy formation based on the Millennium II N-body simulation into 
the semi-numerical modeling of reionization. The semi-analytical model is used to predict the 
production of ionizing photons, then we use the semi-numerical method to model the reionization 
process. Such an approach allows more detailed modeling of the reionization, and also connects 
observations of galaxies at low and high redshifts to the reionization history. The galaxy formation 
model we use was designed to match the low-z observations, and it also fits the high redshift 
luminosity function reasonably well, but its prediction on the star formation falls below the 
observed value, and we find that it also underpredicts the stellar ionizing photon production 
rate, hence the reionization can not be completed at z ^ 6 without taking into account some 
other potential sources of ionization photons. We also considered simple modifications of the 
model with more top heavy initial mass functions (IMF), with which the reionization can occur 
at earlier epochs. The incorporation of the semi-analytical model may also affect the topology 
of the HI regions during the EoR, and the neutral regions produced by our simulations with the 
semi-analytical model appeared less poriferous than the simple halo-based models. 

Subject headings: Cosmology: Reionization, Galaxy: Formation 



1. INTRODUCTION 

The reionization of the hydrogen gas in the 
universe is a topic of forefront research in recent 
years. The measurement of the cosmic microwave 
background (CMB) polarization by the Wilkinson 
Microwave Anisotropy Probe (WMAP)indicates 
that the reionization occurred at Zreion = 10.6 ± 
1.2, assuming an instant reionization model 



([Larson et al.l 120111 ). We also anticipate more 



precise measurement from the Planck satellite 
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(jMitra et al.l l201ll: lAhn et al.l |2012| ). On the 
other hand, observations of high redshift quasars 
show the presence of Gunn-Peterson trough at 
z ^ 6, which mark ed the end of the hydrogen 



reion ization process (jBecker et al.ll2001t iFan et al 



20021 ). While we still do not have detailed knowl- 
edge a bout the nature of the ionizing photon 
sources. fRauch et"aL ( 1997 ) found that for the ob- 
served luminosity function, high redshift quasars 
failed to produce enough ionizing photons to 
keep the universe ionized before z ~ 4, thus 
the stars probably contributed the majority of 
ionizing photons, but the fractions of different 
stellar populations and the quasars are presently 
unknown. However, as the capability of space 
and ground based optical/infrared telescopes im- 
proved, we are learning more and more about 
the galaxies at the epoch of reion i zation (EoR) 



( Ota et al.ll2010l:IWilkins et al llioiol: iBunker et al 



2OIOI : iLorenzoni et al.l I2OIII lYan et al.l l20n' 
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Bouwens et"aLll2012l: iBradlev et aLll2012l ). More- 
over, hundreds of Gamma-Ray Bursts (GRB) 
have been detected by th e SWIFT and Fermi 
spacecraft (jSalvaterra l2012h . many of these are 
at high redshifts, with the highest ones, e.g. GRB 
090429B and GRB 090423 at z > 8, and they pro- 
vide additional informat ion on the star for mation 
history during the EoR ( Ishida et al. 201ll ). 

The 21cm emission from high redshift in- 
tergalactic medium (IGM) could potentially al- 
low us to observe the epoch of reionization di- 
rectly, providing 3-dimensional information about 
the ev olution and morph ology of reionization 



process ljMadau et al.l 119971 ). However, it is dif 



ficult to detect this signal with the presence of 
strong foregrounds. Attempt s have been made 
with the EDGES experiment (iBowman fc Roge^rsl 
20101). and GMR T EoR search (|Paciga et al...20l'ir 
Pen et al. l2009l) . Several low frequency radio tele- 
scope arrays have been or arc being built to detect 
this signal, including the 2 1 CM A El, the Mileur a 
Wide Field Array fMWA'l dJoudaki et al.l l201ll) 



the L ow Frequency Arr ay (LOFAR) 
2010l ). and PAPER (IParsons et al 



Harker et al 
20101) 



In 



the future, HERA (jPurlanetto et al.l l2009l) and 



the S quare Kilometer Array fSKA) (jSantos et al 



20111 ) may provide even more powerful observa- 
tional probes. 

Detailed modeling of reionization is required 
to interpret the various observational data. An 
accurate numerical model must include treat- 
ment of gas dynamics, chemistry, feedback pro- 
cesses, and especially, radiative transfer of the 
ionizing photons. Many groups have conducted 
radiative transfer simulations to study the EoR 
( Gnedinll2000HRazoumov et al.|[2002l: ICiardi et al ' 



20031: ISokasian et al.l l200lt iMaseUi et al.l I2OO3F 



McUcma et al.'2006':'McQ uinn et al.ll2007tlTrac fc CenI 
2007: .Altav et al...200& Aubert fc Tevssiei) I2OO8I: 



Finlator et al.l I2OO9I : IPetkova fc Springell \200^ 
However, such simulations have extremely high 
demands on the computing resources. High res- 
olution is required to resolve low mass galaxies, 
which may contribute a significant fraction of ion- 
izing photons, but as the typical size of the ionized 
regions at the end of EoR is expected to be tens 
of comoving Mpc, a large simulation box is also 
required to statistically sample the distribution of 



HIT regions. This large dynamic range put severe 
demands on the computational cost of the simu- 
lation. Furthermore, since our knowledge about 
high redshift galaxies is still very rudimentary, we 
need to explore a large range of parameter space, 
but the high computational cost make this diffi- 
cult or impossible to do. For these reasons, it is 
very worthwhile to make simpler but faster mod- 
els to gain some physical insights and explore the 
parameter space. 

Inspired by the results of numerical simula- 
tions with radiative transfer computation of ion- 
izing photons, an analytical m odel known as the 
bubble model" was d eveloped (jZaldarriaga et al ' 



2004 iFurlanetto et ah .2004bija|). It uses an ex- 



cursion set formulation of structure formation to 
study the size distribution of ionized regions and 
the induced 21cm emission power spectrum. In 
this model, whether a region is ionized or not is 
determined by the ratio of the number of ionizing 
photons produced locally and in the surrounding 
regions to the number of baryons. This model 
is intuitive, and allows one to calculate the size 
distribution of the HII regions during the EoR 
analytically. This method is then extended to 
the s o called semi-numerical model 
20071: iMesinger fc Furlanettol 12007 



2011 



Zahn et al 



Zahn et al 



Mesinger et al.l 2011 ). First order pertur- 



^http: / /trend. bao.ac.cn/indcx. html 



bation theory is used to produce the halo distri- 
bution function at any given redshift, then the 
star formation rate (SFR) and the number of ion- 
izing photons produced is calculated. Finally, the 
same ionizing photon-to-baryon ratio criterion is 
used to obtain the ioni zation field from the halo 
field. Zahn et al. ( 201l[ ) compared this algorithm 
with a ray-tracing radiative transfer simulation, 
and found that they are in good agreement, even 
though it only cost a tiny fraction of the comput- 
ing time of the full radiative transfer simulation. 
However, these semi-numerical models do not con- 
sider the galaxy formation process in detail. 

The galaxy formation process is complicated. 
Besides the nonlinear evolution of dark matter 
density fluctuations, it also involves the heating 
and cooling of gas, ionization and recombination, 
formation of molecules and chemical enrichment, 
formation of stars and feedback process. Lim- 
ited by the dynamic range of simulations and our 
knowledge about the complex physics of these 
processes, all numerical simulations have to adopt 
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some kind of phenomenological approximation. In 
the so called semi-analytical modeling approach, 
galaxy formation processes, esp. star formation, is 
simulated by following a set of prescriptions based 
on the property of the dark matter halos, and 
after specifying some model parameters, it can 
be used to predict the observational properties of 
galaxies, such as luminosity function, stellar age 
and metal distribution, etc. at different redshifts. 
The model parameters are adjusted to fit the ob- 
servations. Comparing with hydrodynamic sim- 
ulation, the major advantage of semi-analytical 
model is that it has much lower computational 
cost, so a large range of parameter space can be 
explored without repeating the whole simulation. 
Semi-analytical modeling has become a power- 
ful tool of cosmological investigation, and over 
the years the models have been developed to in- 
corporate more physical details, and to provide 
predi c tions of more observati ons! White fc Frenk 
19911 : iKauffmann et all Il993t ICole et all ll994F 



We use the sem i-analytical model developed by 
Guo et al. ( 201l[ ). which was based on the Mil- 



Cole. 



2003 



2007 



lennium and Millennium I I N-body simulations 
(jBovlan-Kolchin et al ] l2009l) . 3y modeling a num- 
ber of physical processes in a plausible way, and af- 
ter tuning the parameters this model successfully 
reproduced many observables, such as the lumi- 
nosity function, star formation rate etc. We then 
use the density field from the same Millennium II 
simulation, with the baryon density tracing the 
the dark matter den s ity, an d the semi-numerical 
model of Zahn et al] ( 201lh to calculate the evo- 
lution of ionization fraction. 

This paper is organized as follows. In section 
2 we briefly review the semi-numerical algorithm, 
in section 3 we introduce the semi-analytic model 
and describe our improvements. In section 4 we 
show our result of the rcionization history and 
morphology. In section 5 we discuss our param- 
eter space and summarize our results. 



Kauffmann et al.lll999 : Somerville fc Primackll999t 2. Semi Numerical simulation 



et al.ll2000l:ISp ringel et al.ll200ll:lHatton et al 



Croton et al.. 2006: Bower et al.l2006[ 



Guo fc Whitdl2009l : IWeinmann et al 



■De Lucia RlaizrS l e semi-numerical mod el of reionization ljZahn et al 
\20l(t . 



Recent semi-analytical models can fit the ob- 
servational data such as the luminosity and mass 
distributions of galaxies and quasars, the history 
of star formation and quasar evolution, and also 
the correlation of a number o f observable prop- 
erties of galaxies and qua sars ( Bower et al. 20061 
De Lucia fc BlaizotI 120071 ). Semi-analytical mod- 
eling has also been u sed to study h i gh re dshift 
universe, for example, iRaicevic et al.l (|201l[ ) used 
the Durham GALFORM model to investigate the 
rcionization history. They compared the ioniz- 
ing photon number with the hydrogen atom num- 
ber in different models of initial stellar mass func- 
tions, and found that the stellar components were 
enough to ionize the universe at z ^ 10. 

In the present paper, we introduce the semi- 
analytical modeling of galaxy formation into the 
semi-numerical model of rcionization. The semi- 
analytical model provides more detailed infor- 
mation about galaxy formation and the pro- 
duction of ionizing photons. More importantly, 
this allows us to compare the EoR models and 
probes with the observations at lower redshift, so 
that eventually a self-consistent model of galaxy 
formation and rcionization could be developed. 



^If^,^^^^^^ ■ ■ I I T 

I2OIII : iMesinger et al.ll201l[) is an extension of th e 
analytical bubble model (jPurlanetto et al.ll2004bf ) . 
It is assumed that before the completion of rcion- 
ization process of the IGM, a number of HII re- 
gions will appear in the IGM, and they are pref- 
erentially located in regions of higher densities, 
because in such regions structures formed earlier. 
The mass of each HII region is proportional to 
the number of ionizing photons produced within 
the region, the criterion of ionization of the region 
is fcoU > C~^i where ^ is the efficiency parame- 
ter, fcou is the fraction of mass in collapse object. 
It could be written as C = fescf*N^/{l + Nrec), 
where fesc is the fraction of ionizing photon which 
escaped the halo into the IGM, is the fraction 
of baryons in stars in the halo, is the number 
of photons emitted per baryon in the stars, and 
Nrec is the average number of recombinations. We 
can rewrite the rcionization criteria as 



(1) 



where 



4(m, z)^5e- ^/2K{C)[cJLn - o\m)Y^I'' (2) 

and i^(C) = erf^^(l — Q. 6m is density fluctua- 
tion and 5x is the rcionization criteria. A point 
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in a region of mass m is marked as ionized if 
and only if the scale m is the largest seale for 
which the condition Eq. ([ij is satisfied. The mass 
function of the HII regions can then be obtained. 
For normal stars, it was estimated that a plausi- 
ble set of parameters are /* ~ 0.2, ^ 3200, 
fesc ^ 0.05, Nrec ^ 2, SO ( ^ 10. However, there 
are large uncertainties in all of these parameter 
values. For example, most low redshift observa- 
tions indicate fesc < 5%, but the escape fraction 
at high redshifts might be much larger. Obser- 
vations of galaxies at z = 1 — 3 indicate a broad 
range of escape fr action values from a few percent 
to tens of percent ( Shapley et al. 2006t Siana et all 



2007t iGrazian et al.ll201l[ llwata et al.ll2009[) . 71s 
also very uncertain, and at pres ent vastly differ- 
ent ch oices of ( value is possible ( Furlanetto et al.l 

2004bl). 

Zahn et al.l (|201l[) and iMesinger et~all (|201l[ ) 



generalized the analytical bubble model to what 
they called semi-numerical simulations, and the 
latter are publicly avaiable known as the 21cmFAST. 
The procedure for such calculation is: 

1. Creating the linear density and velocity 
fields; 

2. find halos from the density field; 

3. reallocate halo position by first order 
perturbation theory; 

4. generate ionizing field by the equation: 

C * rrigai > ruH 

In step 2 and 4, the formation criteria are 
checked from large scales down to each single cell 
of the simulation box to flag a halo or an HII re- 
gion. Once the criterion is satisfied, a halo is gen- 
erated (in step 2) or an ionized region is marked. 
For the ionization, one could either flag all pixels 
inside the region as ionized, or only flag the cen- 
ter pixel as ionized. Obviously, the latter is much 
faster, but the results tur ned out to be simila r with 
no significant difference. I Zahn et al.l (j201ll ) com- 
pared their result with a radiative transfer simu- 
lation with the same initial conditions, and found 
the result of the semi-numerical simulation was a 
good approximation of the radiative transfer sim- 
ulation. Thus, the semi-numerical algorithm cap- 
tures the bubble topology and rcionization history 
quite well with moderate amount of computation. 



In this paper, we further improve the semi- 
numerical model by implementing a more detailed 
model of ionizing photon production based on the 
semi-analytic model of galaxy formation. Both 
the semi-numerical model of reionization and the 
semi-analytical model of galaxy formation are effi- 
cient in computation, so our model still allows rel- 
atively quick exploration of the parameter space. 
Moreover, this approach also allows us to investi- 
gate how the physical processes affect reionization 
history, and to constrain the reionization model 
parameters with galaxy observations. 

To include the physical process in the galaxies, 
we rewrite the ionizing criteria as 



^■y fesc 
1 -\- JS^j'ec 



> NiGM, 



(3) 



where Njgm is the hydrogen number density 
in the IGM. The Millennium II simulation is a 
pure dark matter simulation, we assume that the 
baryon density traces the dark matter density, 
Pb = Pd* ^bl^m- For the rcionization simulation, 
we smooth the density field on to 256'^ grid. We 
calculate for each galaxy by relating its UV 
luminosity with the star formation rate, and by 
integrating this luminosity through its formation 
history we could get the total number of ionizing 
photons. Our algorithm is: 

1. convert dark matter density field in 
Millennium II simulation to IGM density 
field; 

2. located galaxies from semi-analytical model 
into simulation box and calculate ionizing 
photon number; 

3. generate the ionizing field by Eq.Q. 

In step 2, we assume that for ordinary Pop I 
stars the total number of ionizing photon s is sim i- 
lar to the value given in iFurlanetto et al.l (|2004bf ). 
However, we also tested how the stellar population 
affects the ionizing history. In our model C (see 
next section) we treat the luminosity as a func- 
tion of metallicity. I n step 3, for e ach p ixel we use 
the same method as Zahn et al] ( 201l[ ). checking 
the ionizing criteria from a large scale compara- 
ble with the simulation box and step downwards. 
Once the criteria is satisfied we flag all pixels in- 
side the region as ionized, or if not, we move into 
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a smaller radius to repeat the above process, till 
the region is ionized or we have reached a single 
pixel. 

The effect of recombination at high redshift is 
quite complicated, since the purpose of this pa- 
per is to develop a relatively fast method to study 
the ionizing history and the morphology of HII 
regions, here we shall make a simple assumption 
of a constant num ber of recombinations here (see 
Yue &: Chenll2012l for an example of sophisticated 
modeling of the evolution of recombination rate) . 
Here we adopt N^ec = 2. 

3. Semi- Analytical Model 

In the hierarchical structure formation sce- 
nario, the dark matter halos grow by accretion 
and merger, and galaxies form within dark mat- 
ter halos by the radiative cooling of the bary- 
onic gas. In semi-analytical models, one fol- 
lows the evolution of each dark matter halo, and 
apply a set of rules to describe the gas cool- 
ing, star formation and feedback for each halo 
without actual simulations. The properties of 
the galaxies within each halo, such as the stel- 
lar mass, the age of stellar population, distribu- 
tion within halos, and the amount of cold and 
hot gas and metal abundance can be followed. 
The dark matter halo merger tree is generated ei- 
ther by Monte-Carlo simulations( White fc Frenk 
ll99ll : lKauffmann et al.lll99i ICole et al.lll994h . or 
by using an N-body simulation. With improve- 
ment in computing power, in most recent works 
the latter approach is adopted ( ^Bowcr et al."200 6t 
2006tlDe Lucia fc Blaizot.2007uGuo fc 



ficient mass resolution that is currently available. 
It assumed a KCDM cosmology with parame- 
ters based on combined analysis of the 2dFGRS 



(IColless et al.l 12001 ) and the first year WMAP 



20031 ). The parameters are: 
0.045, VLa = 0.75, n = 1, 



data (jSpergel et al 
fl™ = 0.25, 1]& = 
(Ts = 0.9 and = 73kms^^Mpc~^. These pa- 
rameters differs slightly from the more recent best 
fit values (jLarson et al.l 120111 ) but the relatively 
small off-set are not significant for most of the 
issues discussed in this paper, as there are much 
larger uncertainties in star formation and reioniza- 
tion parameters. There are 2160^ particles in MS 
II, the box size is 100Mpc//i, the softening length 
lkpc//i, and the particle mass is 6.9 * 10^Mq//i. 

IGuo et aL ( 2011 ) developed a semi- analytical 
model based on the MS II. This model is an exten- 
sion and improvement of ear lier models based on 
the Millennium sirnulation ( Springel et al. _ 2005t 
Croton et all 120061: lOe Lucia fc Blaizotll2007t )"~ln 



Croton et al 



20091 : IWeinmann et al.ll2010l ). 

For galaxy formation modeling during the 
epoch of reionization, high resolution simulation is 
required. Indeed, it is believed that the first stars 
form in halos of lO^M©, and the first galaxies 



Whit 



this model, central galaxies, satellite galaxies in 
subhalos and orphan galaxies which had lost their 
subhalos are distinguished. The baryonic content 
of the galaxies is split into five components, in- 
cluding a stellar bulge, a stellar disk, a gas disk, a 
hot gas halo, and an ejecta reservoir. In addition, 
intra-cluster light is also included. The model keep 
track of various processes involved in galaxy for- 
mation, including gas heating and cooling, evolu- 
tion of the stellar and gas disks, star formation, 
supernova feedback, gas striping in groups and 
clusters, merger and tidal disruption, bulge forma- 
tion, black hole growth and AGN feedback, metal 
enrichment, and photo- heating of the pregalactic 
%as by the UV background afte r reionization. A 



Chabrier IMF (jChabrieij 120031 ) is adopted, this 
IMF has fewer low-mass stars than the Salpeter 
IMF. The mo del use the stellar popula tion syn- 
thesis model of lBruzual fc Charloti ( 20031) . and the 



have masses about 10 Mq (c.f. iBarkana fc Loeb 
20011 ). At the same time, the model should also 
contain enough volume such that a large num- 
ber of neutral or ionized regions can be included 
in the simulation at the EoR. In this work, we 
use a model b ased on the Millennium II simula- 
tion (MS II) (jBovlan-Kolchin et al.|[2009l ). This 
is an e xtension of the earl ier Millennium simu- 
lation (jSpringel et al.l |2005[) . and is among the 
largest cosmological N-body simulations with suf- 



dust e xtinction model developed in lGuo fc White 
(|2009f ). It predicts observable properties such as 
the stellar mass function, mass-size relation, dis- 
tribution of galaxies within cluster and group ha- 
los, morphological type, black hole mass - bulge 
mass relation, low redshift galaxy luminosity func- 
tion for different observing bands, stellar mass- 
halo mass relation, cold gas metallicity, galaxy 
color distribution, Tully-Fisher relation, satellite 
luminosity function, autocorrelation function for 
different masses, etc. Some of these predictions are 
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compared with observations to fix the model pa- 
ramete rs. Compared with earhe r semi-analytical 
model (jPe Lucia fc Blaizotll2007l ). which overpre- 
dicts the abundance of galaxies with mass near 
or below IO^Mq when applied to the MS II, this 
model improved the treatment of a number of 
physical processes, including the treatment of su- 
pernova feedback, reincorporation of ejected gas, 
sizes of galaxies, treatment of satellite galaxies 
which are outside the i?2oo but belongs to the 
same group, and environmental effects. The model 
is adjusted to best fit the observational data on 
galaxies at 2: ~ 0, as there are much more high 
quality observational data available at low red- 
shifts, and they are subject to less selection effect. 
In this paper we adopt this model as our model A 
and apply it to the semi-numerical simulation of 
reionization. 

While the predictions of this semi-analytical 
model are in good agreement with many obser- 
vations at z < 1, and the abundance of luminous 
galaxies are also consistent with observations at 
z < 3, the predicted high rcdshift star formation 
rates are systematic ally lower than t he observed 
values (see Fig. 22 of IGuo et al. 2011 and discus- 
sions there). This is not simply a problem with 
this particular model, for it is well known that 
if the observed star formation rate is integrated 
over redshifts, the luminosit y function of galaxie s 
would b e overpredi c ted (c .f. IWilkins et al. 2008} ). 
In fact, iGuo et al. 1 (I2OIII) argued that even the 
present model might have overpredicted the abun- 
dance of low mass galaxies at high redshifts. How- 
ever, as we shall see in the next section, with this 
model the ionizing photon production rate might 
be too low to allow the universe to be reionizcd at 
z>6. 

There have been other semi-analytical studies 
which tries to accommodate more of the high red- 
shift SFR observations. For example, in most 
semi-analytical model it is assumed that major 
mergers can trigger starbursts, so that most gas 
in the merging galaxies collap s e and form stars in 



a short time. iRaicevic et al.l (j201lh argued that 



in such a scenario the stars are much more heavy 
than stars formed continuously. So they assumed 
a top heavy initial mass function for the stars 
formed in major mergers, and the number of ioniz- 
ing photons would be increased dramatically. Ma- 
jor merger are more frequent at high redshifts, so 



this model can be expected to increase the pro- 
duction of ionizing photons by a factor of 5 to 10 
during the EoR. 

To achieve a higher photon production rate dur- 
ing the EoR, we may c onsid er a more top heavy 
IMF as Raicevic et al.l ( 2011 ) did. Strictly speak- 
ing, once a semi-analytical model is set, its model 
assumptions or parameters such as the IMF should 
not be altered, because the feedback would change 
subsequent star formation, and if one runs the full 
semi-analytical model with the new assumptions 
and compare the results with observations, a dif- 
ferent set of parameters would be obtained. How- 
ever, as we noted above, at present there is a con- 
flict between the low redshift abundance and high 
rcdshift star formation rate, which can not be eas- 
ily solved. So to illustrate the effect of change, 
here we will just make a simple modifications of 
the IMF, and examine its impact without consid- 
ering the feedback. I n our model B, we ma ke an 
assumption similar to iRaicevic et al.l (|2011l ). that 
major merger will trigger star burst, the IMF in 
star burst galaxies are very top heavy with an in- 
dex of 0, mass range of (0.15 — I2OM0). By doing 
this, more than half of the newly formed stars are 
massive ones, and the number of ionizing photon 
is about 10 times than model A. 

As the model B may be too extreme, we also 
consider a model C. As noted above, during the 
epoch of reionization, many of the stars arc formed 
in metal free or very low metallicity regions, which 
may have a top heavy IMF. Furthermore, even 
for the same mass, the metal poor stars produce 
more ionizing photons. Here, we do not assume 
the flat IMF as in model B, but instead consider 
a mor e moderate model as suggested in ISchaerei 
(|2003l ). who examined the spectral properties of 
the ionizing continua, the Lyman break and re- 
combination lines in star bursts and constant star 
formation phase, for metallicity from up to solar 
metallicity. We adopt a Salpeter IMF for all galax- 
ies, attributing very massive stars in the range of 
(1 — 5OOM0) with index of 2.3, and mass range 
of (1 - lOOM©) with index of 2.55 for the galax- 
ies with higher metallicity. Applying these results, 
the procedure of our simulation with model C is 

1. The ionizing luminosity for a galaxy at a 
given redshift is determined by its 
metallicity, which is determined by the 
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metallicity of the cold gas at last redshift. 

We interpolate t he property of t he stellar 
population from ISchaered (|2003f ) and 
Bruzual fc Chariot (,200sl l. to obtain the 



luminosity; 

3. Integrate the luminosity in the past to get 
the number of ionizing photons; 

4. Generate ionizing field with the number of 
ionizing photons from step 3 and IGM 
density. 

We may estimate the mean metallicity of the 
newly formed stars from the the metallicity in the 
cold gas of the halo. However, in the original semi- 
analytical model, newly formed halos (i.e. without 
progenitors), regardless of redshift, will always be 
metal free, as the metals are not transferred out 
to the IGM in that model. In reality, the IGM will 
be polluted by galaxy winds, which blow the en- 
riched gas out of the halo, and Lyman alpha forest 
observations show that the IGM has alrea dy been 
contaminated at z > 3 ( Schave et al.ll2003l ). Thus, 
even the newly formed halos should contain some 
amount of metals. As a remedy of this, we set 
the metallicity in the cold gas of the newly formed 
halos as the average value of the whole simula- 
tion box, which is calculated by the metal pro- 
duction in the original semi-analytical model. In 
Fig. [1] we plot the star formation rate of different 
metallicities according to this model, the Pop I, 
II, and III stars are represented by the blue dot, 
green dash and red dot-dash curves, respectively, 
and the black solid curve represents the total SFR. 
While initially (at very high redshift ) the Pop III 
had the largest SFR, this was soon surpassed by 
Pop II {z ~ 16), which was in turn surpassed by 
Pop I at a fairly high redshift {z ^ 12). The for- 
mation of Pop III stars declined at z ~ 10, and 
fluctuates a little at low redshift but stays at low 
value. In this model we have assumed the metals 
are uniformly distributed, the non- uniform distri- 
bution of metals may allow a higher production 
rate of metal free or metal poor stars. 

In Fig. [2J we plot the evolution of the ionizing 
photon cumulative production rate due to stars of 
different metallicities in our model C. We see that 
during most of EoR, the contribution from the Pop 
II stars dominates. The Pop III stars could make 
up major contribution only at fairly high redshifts. 




Fig. 1. — The SFR of different metallicities in our 
model C. Here the criteria of Pop II and Pop HI 
stars are -^Zq and 10"'' Z© respectively. 

but levels off at z ~ 10 in this model. The Pop I 
stars have a higher SFR, but do not make up the 
major contribution to the ionizing photons. 

4. Results 

Let us first take a look at the global production 
rate of ionizing photons in our model. In Fig|31 
we plot the ratio between the number of ionizing 
photons from stars and the number of hydrogen 
atoms per comoving volume, here the escape frac- 
tion and recombination number are not included. 



Thus, if we assume fe 



0.15 and iV^, 



2, the 



ionizing photo-to-baryon ratio has to be greater 
than 20 for the Universe to be reionized. The blue 
dash line shows the result for assuming that each 
hydrogen atom in the stars produces on average 
3200 photons during the life time of the star, which 
corresponds to the normal stellar population given 
by our model A. The red dot dash line shows the 
result for our model B, which produce much more 
photons as we assumed a very top heavy IMF 
for starburst galaxies. The black solid line cor- 
responds to our model C, which assumed a mix- 
ture of Pop I, Pop II and Pop HI stars. From the 
figure, we can see that for model A, it is very dif- 
ficult to ionize the universe by z ^ 6, even if large 
escape fractions and low recombination rates are 
assumed. Model B could ionize the universe at 
relatively early time, perhaps z ~ 8. In model C, 
reionization occurred at z ~ 7, which is still much 
later than the z,.e = 9 given by the WMAP data 
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Fig. 2. — The number of ionizing photons of differ- 
ent metalhcities in our model C. Here the criteria 
of Pop II and Pop III stars are -^Zq and 10~'^Zq 
respectively. 



( Larson et al.ll201l[ ) 



We may also compare the galaxy UV luminos- 
ity function with observations at high rcdshifts. 
Figs. |3][5] show the 1500A luminosity functions 
predicted by our model A and model C at z = 7.88 
and z = 7.27 respectively, and for comparison 
we also plot the measured luminos ity function at 
I6OOA from lBouwens et all (|201lh . We see that 
our model A is in excellent agreement with the lu- 
minosity function, even though its prediction on 
SFR falls below many observations. Our model C 
overpredicts the luminosity function slightly, but 
given the large uncertainties in the data at present, 
the deviation is still acceptable. 

With the semi-analytical model of ionizing pho- 
ton production in place, we can now use the pro- 
cedure described in Sec 15] to simulate the reioniza- 
tion process. Figure |6] shows the reionization pro- 
cess of our simulation box for our model C. The 
black colored regions are neutral, while the white 
ones are ionized. The rcdshifts are from the upper 
left to bottom right 14.9, 11.9, 10.9, 9.3, 7.88 and 
7.27 respectively. We can see that the universe is 
mostly neutral before redshift 12, but with some 
bubbles of HII region. Around redshift 9, the bub- 
bles began to overlap, and by z = 7, most of the 
universe is ionized. 

How does this semi-numerical model compare 
with a model without the semi-analytical model- 
ing? This depends somewhat on the particular 



Fig. 3. — The ratio between number of ionizing 
photons from stars and number of hydrogen atoms 
for the three models we considered. 



model we use as well as the redshift. However, 
in some cases the difference can be quite obvi- 
ous. In Fig 17] we plot the ionization at redshift 
z = 6.7 for our model A. We have here chosen 
to show model A for comparison because model A 
uses the original semi-analytical model. As we can 
see from the figure, while the large scale distribu- 
tion of the HI and HII regions are similar, we see 
that in the model without semi-analytical model- 
ing, there are more small HII regions which are 
absent in the model with it. In the model without 
semi-analytical galaxy formation, the number of 
photons are predicted from the dark matter halos, 
as a result, some halos which do not contain stars 
were also assumed to produce ionizing photons. 
This would not affect large HII regions, but in the 
neutral regions, small HII regions are produced 
around these small halos. This can significantly 
affect the topology of the ionization map, giving it 
a more poriferous appearance. In the model with 
semi-analytical modeling, on the other hand, only 
galaxies contribute to ionizing photons, while the 
starless minihalos do not contribute, so the HI and 
HII region are more well separated with a mono- 
lithic appearance. 

In FiglHlwe plot the evolution of the global ion- 
ization fraction for the different models, where we 
assumed fesc = 0.15 and N^ec = 2. These re- 
sults are consistent with the ratio between ion- 
izing number and hydrogen number in figure E] 
For model A, the reionization happened very late. 



8 



- 8, Bouwens(2011) 

- 8, Model A 

- 8, Model C 





Fig. 4. — The luminosity functio n at z = 7 



The data points are from Bouwe ns et al.l (|201l[ ) 



The prediction of our model A and model C are 
shown as blue dash line and black solid line re- 
spectively. 



S. -3 



~1. Bouwens(2011) 

- 7, Model A 

- 7, Model C 




Fig. 6. — Ionizing map of our model C, with pa- 
rameters /esc ~ 0.15 and N^ec — 2. From upper 
left and bottom right, the red shift are 14.9, 11.9, 
10.9, 9.3, 7.88 and 7.27. Black and white region 
corresponding HI and HII region, the pixels could 
only be total neutral or total ionized. 




Fig. 5. — The same as Fig. lU but for z = 1.21. 

even at z = 6, the universe is still not ionized in 
this model. In model B, the ionization occurred 
at z '--^ 8. In model C, the universe is fully ionized 
at z '--^ 7, and it reach ]hii = 0.5 at z ^ 9. 

5. Summary and Discussion 

In this paper we incorporated a Millennium II 
based semi-analy tical model of galaxy formation 
( Guo et aI1l2011 ) in the semi- numerical simulation 
of reionization, which was shown to be a good ap- 
proximation to the radiative transfer computation, 
but with much less computational cost. Star for- 
mation rate is given by the semi-analytical model, 
it is then used to compute the production rate 
of ionizing photons in each galaxy. As this semi- 



Fig. 7. — Comparison of ionization map produced 
by the semi-numerical simulation without (left) 
and with (right) semi-analytical models. 



analytical model is exclusively designed to best fit 
the observations at z '--^ 0, there is some discrep- 
ancies between its predictions and high redshift 
observations. This model, denoted our model A, 
predicts an SFR lower than the observed values at 
high redshifts, but it does reproduce the observed 
UV luminosity function at z ^ 7 — 8, thus clearly 
revealing the conflict between the observed lumi- 
nosity function and SFR at high redshifts. How- 
ever, the ionizing photon production rate for this 
model is too low, such that the Universe can not 
be reionized at z > 6 for an escape fraction of 
/esc — 0.15 and the mean number of recombina- 
tion Nrer ~ 2. 
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6 7 8 9 10 11 12 13 14 15 



Fig. 8. — Evolution of the ionization fraction for 
the different models, here we accept /esc — 0.15 
and iVrec = 2. The vertical green and yehow 
lines mark the reionization redshifts constrained 
by Gunn-Peterson trough of quasars and WMAP 
observations. 

To remedy this problem, we also considered two 
simple revised models. These models assume the 
same star formation rate as the original model, 
but adopted different IMF for star formation. In 
model B, a flat IMF is assumed, this model can 
produce about 10 times of the ionizing photons of 
model A, and the universe is reionized at z = 8. In 
model C, we assumed a top heavy IMF but with 
similar slopes, and we also consider the production 
rate of ionizing photons as a function of mctallicity 
of the newly formed stars. The metallicity of the 
newly formed stars is determined by the metallic- 
ity of the cold gas in the halo, which was given in 
the original semi-analytical model. We also test 
other definition of metallicity, for example, the ra- 
tio between the increase of the total mass of metals 
in stars and the total mass of new formed stars, 
i.e. the averaged metalicity in new formed stars. 
This definition reduces the total amount of the 
photons by 10% - 20% but won't affect our main 
conclusions. And we assumed a metallicity floor 
in the IGM as given by the average metallicity of 
the simulation box. This model predicts that Pop 
II stars are the major sources of ionizing photons 
during the EoR, even though the Pop I SFR is 
higher. The predicted UV luminosity function at 
z ~ 7—8 of model C is only slightly higher than the 
observed value. The universe is ionized at z ~ 7 
in this model. 



Armed with the semi-analytical model, we sim- 
ulate reionization history with the semi-numerical 
method. The addition of the semi-analytical 
model could significantly affect the topology of the 
neutral and ionized regions. While it is somewhat 
dependent on the particular model and redshift, in 
some cases the difference is quite large. Without 
the semi-analytical model of galaxy formation, we 
may incorrectly assume ionizing photons were pro- 
duced by starless minihalos in HI regions, and this 
may produce a more poriferous appearance. With 
the semi-analytical model, the ionizing photons 
come only from galaxies, not from all minihalos, 
so the Hll regions appeared more monolithic and 
well separates from HI regions. 

There are some simplifications wc made in this 
research. For example, we have adopted simple 
constant mean value for escape fraction and num- 
ber of recombinations. In reality, both of these 
may evolve with redshift, or even vary with differ- 
ent galaxies. We have assumed a simple two phase 
model of IGM, in any given point it is either fully 
neutral or fully ionized, partial ionization is not 
consider, nor do we consider the effect of quasar 
formation. However, within the plausible range of 
the parameters, these effects are unlikely to change 
our result qualitatively. As the aim of this paper 
is to illustrate the use of semi-analytical model of 
galaxy formation in semi-numerical modeling, we 
do not try to model these complications. 

There arc certain limitations in our present 
work. Perhaps the most important one is that our 
model B and C are not completely self-consistent. 
We have assumed that they have the same star 
formation rate evolution as given by model A, 
but with modified IMF and ionizing photon pro- 
duction. However, once the IMF of the model is 
changed, the star formation history would not be 
the same, because the feedback effect would be 
different, and that would change the subsequent 
accretion and star formation. If wc hope to have 
a similar star formation history with this mod- 
ified IMF, we must assume a lower strength of 
feedback. Furthermore, the reionization itself may 
also affects the galaxy formation process. To be 
really self-consistent, one has to modify the semi- 
analytical model itself. This will be the next step 
of our research. 

As wc discussed above, the semi-numerical sim- 
ulation combined with semi-analytical galaxy for- 
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mation model provide a good approximate of 
reionization history. It allows us to investigate 
large parameter ranges in short time, and more 
importantly, connects the observations on galax- 
ies to the reionizations. This approach could help 
us better understand the galaxy formation pro- 
cess. 
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